Optimal approximation of harmonic growth clusters by orthogonal polynomials 
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Interface dynamics in two-dimensional systems with a maximal number of conservation laws gives 
an accurate theoretical model for many physical processes, from the hydrodynamics of immiscible, 
viscous flows (zero surface-tension limit of Hele-Shaw flows, pQ), to the granular dynamics of hard 
spheres 2 , and even diffusion-limited aggregation 3 . Although a complete solution for the con- 
tinuum case exists [4, 5 , efficient approximations of the boundary evolution are very useful due 
to their practical applications [6 . In this article, the approximation scheme based on orthogonal 
polynomials with a deformed Gaussian kernel [7] is discussed, as well as relations to potential theory. 

PACS numbers: 05.70.-a, 02.10.Yn, 02.30.Tb, 05.30-d, 05.40.-a, 05.50.+q 



Introduction - A large class of two-dimensional pro- 
cesses, both deterministic and stochastic, have been 
mapped to a powerful mathematical model known as har- 
monic (or Laplacian) growth [ 8| t9| fTU l fTT] fT2 | fT3 l fT4| fT5 | 
[T6] . In this model, a domain D+(to) grows in time t, i.e. 
its normalized area to = 7r _1 J d dxdy (areas in units of 
7r throughout) has a linear dependence to = Qt, with Q 
the pumping rate. 

This prescription is obtained as a consequence of the 
complete formulation as an exterior problem: if T rep- 
resents the boundary of D+ and D_ its complement, we 
seek the harmonic field p (which may represent pressure, 
concentration, etc) and the velocity field v, such that: 



-d n p, P = on T, 
-Vp, Ap = onD_, 
- log |z|, z — > oo G D- 



(i) 



with v n , d n p their normal components on T, respectively. 




FIG. 1: Simply connected droplets and equipotential lines for 
the interior and the exterior domains. 



An equivalent formulation of the interior problem was 
studied by S. Richardson [17], where the roles of exterior 
and interior domains are interchanged: p is harmonic 
inside up to a finite number of isolated logarithmic sin- 
gularities (sources and sinks). Assume for simplicity a 
single point source at the origin z = with pumping rate 
Q, or, in terms of the pressure p, Ap = & S 2 (z). Richard- 
son proved that for any function h G L 1 (.D + ) harmonic 



in a neighborhood of D+ we have 



d(h). 
dt 



Qh(0), 



(2) 



where (ft) + = f D , f , h(z)dxdy. Therefore, taking h(z) 
z k we have that the interior harmonic moments 



Vk 



* J D + (t) 



z dxdy, 



k=l,2,... 



(3) 



are preserved while the area grows linearly with t. 

Mapping the interior problem by means of a simple 
coordinate change z = 1/z and changing the source to a 
sink (Q — > —Q) implies that in the original exterior prob- 
lem for D- , G D- , the exterior harmonic moments 



tk 



1 

kir 



dx dy. 



1,2,... (4) 



D-(t) 



(the interior moments of D_ now) are preserved by the 
evolution and dto/dt = Q [18], and we set Q = 1. 

For a simply connected droplet D + , the solution of the 
problem is easily expressed through the conformal map 
taking the exterior of the unit disk \(\ > 1 univalently 
onto D-, and matching the points at infinity, 



z = f(C)=rC + J2 u kC k , r>0, 



(5) 



fc=0 



with respect to which the pressure and normal compo- 
nent of boundary velocity read 



p(z) = -log\f- 1 (z)\, v n = \d z f- 1 (z) 



(6) 



As shown in [5], conservation of exterior harmonic 
moments is equivalent to the statement that there is a 
canonical transformation from the variables z.z^, 



* = /(0, «* = /(r 1 ), 



(7) 



to the variables to,logC, i.e. C 

note the alternative formulation dz A dz^ 



dz dz* 
dC dt 



z$ dz 
dC dt 



= 1. We 

d log (" A dto . 



In this paper, we describe an optimal approximation 
procedure for the boundary T(i), as well as for the 
Cauchy transform of D + , 



CdM= 1 -I 

k Jd + 



dx'dy' 

z — z' 



iy\ zeD_. (8) 



An example - A superficial analysis of the model de- 
scribed here would lead to the conclusion that, given the 
integrability of the system, specifying a set of exterior 
harmonic moments {t^} at a given, initial area Tq, should 
yield a solution for to > T without major analytical dif- 
ficulties. In reality, reconstructing the conformal map ([5} 
which corresponds to the data To, {t^} is typically a very 
difficult problem. There are few known correpondences, 
between classes of moments {tk} and classes of shapes (pj 
(see [7] for examples). An illustrative case occurs when 
the exterior harmonic moments form a simple sequence 



tk 



P 



fe = l,2, 



(9) 



where (3 G M + and a G C \ {0} are given parameters. 
Defining two radii 



Ri 



P, 



Ri 



1 + 2/3 



(10) 



we have two completely different cases: if \a\ + R\ < 
R2 then we have a doubly connected domain bounded 
by circles and for \a\ + R\ > R2 the domain is simply 
connected given by an exterior conformal map of the form 



f(()=r( + u 



C-A' 



r > 0, \A\ < 1, 



(11) 



(see FIG. p}. The parameters of /(£) are related to the 
deformation parameters as follows: setting u = v/A, 



P 



to- 
r 
A 



TV 

I 2 ' 



where to 



-*2. 



V 

A 

|2\2 



vA 



(12) 



1-WI 2 ' 



v\ 2 /(l - \A\ 2 ) 2 is the area. If \a\ + R 1 > 



R2 then the above equations have a unique solution for 
r, v and A in terms of (3 and a. 



cSC 



FIG. 2: Laplacian growth for the exterior moments (19} 



Domain approximation via orthogonal polynomials- 



The fact that, even for a rather simple shape (11) the 



correspondence (12) is quite intricate shows the need for, 



and practical value of, efficient approximation methods. 
To that end, we define the following family of orthog- 
onal polynomials: consider the domain specified by the 
exterior harmonic moments {£&} and t$. The function 
defined in a neighborhood of the origin by 



V(z) = J2tkZ k , 



0, 



(13) 



fc>l 



is preserved by the harmonic growth. Now consider the 
function W(z) = \z\ 2 — 25RV(z), which we label confining 
potential, and suppose that 

I \ z \n e ~NW{z) d 2 z <OQ n = 0, 1, 2, . . . 

for all for values of the scaling parameter N > 0. For 
fixed TV, the orthogonal polynomials {Pn (z)} of the 
weight function e ~ NW ( z ) are defined by 

/ Pi N \z)W\z)e- NW ^d 2 z = S nm . (14) 

Jc 

The approximation method presented in this work is 
based on the following statement: as 



00, N — > 00, 



n 

N 



to, 



(15) 



the weighted polynomials \Pn (z)\ 2 e NW ( Z ) (which we 
denote by pn (z) in the following) converge to the con- 
formal measure of the domain D+(i), with support T(t): 



,W(^ = 



(z) = \d z f- 1 (z)\ 



O 



N 



, n,7V- 



00. 



(16) 



The proof of this result appeared first in [7j. We do 
not repeat the entire argument, as it would require too 
much space, but recollect the main ideas: starting from 
the differential equations satisfied by the weighted func- 
tions VJ n (z) = P n (z)e~ NW ( z ^ 2 , with respect to vari- 
ables n/N and z, we integrate perturbatively in powers 
of A/" -1 , and obtain the expression ([7], equation (76)): 



^) - v / [FWexp 



TV 



3? / S(()d( 



where the Schwarz function S(z) is defined by the iden- 
tity S(z) = z, z G T. It is known to have the expansion 



S(z) = V'{z) 



1 r dx'dy' 

7T J D+ Z~Z f ' 



OO 



(17) 



Since the exponent in the asymptotic expression of ip(z) 
vanishes on the boundary T and gives a Gaussian de- 
cay away from it, the weighted polynomials p n (z) are 



described, in the n, N — > oo limit, by the conformal mea- sensitive to the accuracy of the Gram matrix and thus 

requires very precise computation of {g\j} • Then the 



sure (16). However, this asymptotic result says very little 



about the behavior of p n \z) for finite values of n, N. 

In the remainder of this paper, we present numerical 
evidence for the convergence properties of p n (z) at fi- 
nite values of their order. We show that the agreement 
between p n \z) and |[/ _1 (^)] / | is excellent for values 
of n as little as n = 20, and present potential appli- 
cations of this property. One obvious consequence is 
related to reconstruction algorithms of domains in this 
class: assume that the Schwarz function of the domain 
D+ has a branch cut on 7 G D+, with real, positive jump 
function p s (z). Obviously, this may include the case of 
meromorphic functions with poles in z p G .D+, for which 
p s (z) = J2 a p S(z — z p ). Then from the asymptotic result 



llog|pW(z)|->W(s) + ±jf ) 



log\z-z'\ 2 dx'dy', 



(18) 

it follows immediately that the asymptotic distribution 
of zeros of the orthogonal polynomials converges to p s (z): 



1 



N 



hm — > log \z- 

2=1 



f p s (z')log\z-z'\ 2 dz'. (19) 



In a forthcoming publication [19] , we will give detailed 



proofs of (16, 19). In this Letter, we provide a thorough 



numerical analysis supporting the asymptotic results. 




FIG. 3: Localization for the density for n 
formal measure. 



20 and the con- 



density p n (z) is obtained from the polynomial P n (z). 
Of course, the usefulness of this approximation scheme 



relies on the rapidity of the convergence in (16), which 



may not seem to be very promising. However, our nu- 
merical experiment (FIG. [3]) shows that in the example 
(Jol) above the "shape" of the conformal measure (the blue 
curve) is recovered very accurately by the weighted poly- 
nomial density of a degree as low as n = 20. 





FIG. 4: Density plot and contour plot of the localized density 



To illustrate the speed of convergence in (16), the 



Kullback-Leibler divergence or relative entropy of p n (z) 
with respect to the harmonic measure, given by 



hC^wh 



logpi N \f(e ie )))de, (21) 



2Wo V"°l/'(e*)l 
is calculated and plotted on the figure below. 



J N )f 



Simulations and numerical study - Let to and the ex- 
terior harmonic moments be given through the potential 
V(z). To fix the scaling limit, let N(n) = n/t$. For fixed 
n, we have to calculate the entries of the Gram matrix 



9 



(n) 



z i z j e- n/toW{z) d 2 z z, j = 0, . . . , n. (20) 



For potentials W(z) that are converging rapidly enough 
to infinity as \z\ — > 00, the exponentially decaying weight 
makes the planar numerical integration a feasible task. 
The stabilized Gram-Schmidt Algorithm provides the or- 
thogonal polynomials P U5 at(z), which is known to be very 



FIG. 5: The Kullback-Leibler divergence of pn (z) with re- 
spect to the harmonic measure. 



The asymptotic behaviour of the zeroes of orthogonal 
polynomials in the scaling limit (IT5J) was also investigated 
in the particular case ^. Since /(£) is a rational function 
of order two, the Cauchy transform Cjj + (z) (fsl) in the 
exterior domain D_ satisfies a quadratic equation 

A(z)y 2 + B(z)y + C(z) = 0, 

with rational coefficients in z depending on the param- 
eters of /(C)- As an algebraic function, Cd + (z) can be 



analytically continued on a plane with a branch cut con- 
necting up the branchpoints 



£1,2 



A 



At ± 2y/rv 



(22) 



of the inverse mapping f~ 1 (z). This "conjugate elec- 
tric field" created by the uniformly charged domain D+ 
is mimicked by the field generated by the normalized 
counting measure of the zeroes. However, these points 
seem to accumulate along some curve (as opposed to the 
equilibrium configuration in the presence of the back- 
ground potential W(z) - the so-called Fekete points - 
which are distributed asymptotically uniformly). Since 
the asymptotic zero distribution must be real and pos- 
itive, the natural choice is dictated by the Sokhotski- 
Plemelj formula: the critical trajectory 7 is selected by 
the condition that the jump between the two solutions 

y± 



(-B± VB 2 - 4AC)/2A satisfies 

to((y + (z)-y-(z))dz) = 0. 



The critical trajectory can be found by calculating 



$(*) := *R 



(y+( w ) -y-(w))dw 



(23) 



(24) 



and then plotting the contour <&(z) = 0. Three trajec- 
tories are emanating from each branchpoint: there are 
two trajectories that connect z\ and ^2, and the one con- 
tained by the domain attracts the roots. As can be seen 
in FIG. [6j the distribution of zeros (for n = 50) and the 
trajectory 7 are almost indistinguishable. 




FIG. 6: The critical trajectory and the zeroes for n — 50. 



Applications - The method presented in this Letter 
allows to construct optimal approximations with high 
convergence rates for either the boundary or the branch 
cuts characterizing domains from the harmonic growth 
class. This may be used in a number of different situa- 
tions; here we discuss two relevant examples: 

(i) an outstanding problem in viscous two-dimensional 
flows is formation of boundary singularities (cusps). 
They are known to occur for finite values of the normal- 
ized area to, and for many initial conditions [20]. For a 
particular class of such cusps, with local geometry given 



by the scaling x 



y 



2/c+l 



, k = 1,2,..., it is not possi- 



ble to continue the evolution of the boundary beyond the 
cusp formation, and a weaker type of solution is required. 




FIG. 7: Singular shape as a result of Laplacian growth 



The weak solution [21 is based on the equivalence be- 
tween the distribution of zeros of the orthogonal polyno- 
mials and the branch cut of the Schwarz function, (fl9|. 
These two distributions generate the same Newtonian po- 
tential in D- as the uniform distribution on D+ (real 
droplet), so they may be considered as equivalent solu- 
tions before singularity formation, FIG. [7| However, af- 
ter a cusp is formed, smooth (uniform) solutions are not 
possible anymore, while the distribution of zeros of the 
polynomials remains well-defined. The conjecture is that 
this weak formulation will produce solutions which ex- 
plain the famous fingering patterns observed in physical 
realizations of this model. This will be shown in a forth- 
coming publication, where the algorithm will be used to 
construct numerical solutions according to this prescrip- 
tion, and compare with real, physical patterns [22 . 

(ii) another practical application of the results pre- 
sented in this Letter is an efficient algorithm for shape 
(boundary) reconstruction when the domain D+ is given 
through the reduced data to,{tk}. For example, such 
(reduced) representations arise in satellite imaging data 
compression [6]. Shape reconstruction algorithms are 
then needed to find the boundary T, given the set of 
moments to, {£&}, with particular emphasis on good con- 
vergence rates. Since the data to,{tk} is sufficient for 
constructing the family of orthogonal polynomials Pn 
introduced here, we have a boundary approximation al- 
gorithm which gives excellent results already at n = 20. 
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